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1.0  Introduction 


The  subject  of  linear  optimal  approximation  has  received  considerable 
attention  in  recent  years  [4],  [5],  [6],  [7],  [8],  The  subject  of  multi- 
variate approximation  for  scattered  data,  including  optimal  approximations, 
is  reviewed  in  [9].  The  idea  is  appealing  since  the  optimal  approximation 
in  a certain  space  of  functions  minimizes  the  norm  of  the  error  functional 
for  approximations  in  that  space.  When  the  space  is  a Hilbert  space,  the 
computation  of  optimal  approximations  becomes  rather  simple,  in  theory  [2]. 

A known  reproducing  kernel  function  provides  the  representers  of  linear 
functionals  defined  on  the  space.  The  optimal  approximation  satisfies  the 
system  of  equations  obtained  by  requiring  that  the  approximation  be  exact 
for  the  representers  of  the  functionals  being  used  for  the  approximation, 
usually  point  evaluation  functionals. 

In  practice,  it  seems  that  optimal  approximations  have  not  been  used 
very  much.  This  is  perhaps  partly  because  of  a lack  of  experience  with  them, 
as  well  as  the  fact  that  use  of  the  representers  as  a basis  set  for  optimal 
approximation  is  the  analog  of  the  use  of  truncated  power  functions  as  a 
basis  set  for  univariate  spline  approximation. 

The  particular  space  of  functions  of  two  variables  to  be  considered  here 
are  the  Sard  "corner  spaces",  Bp  - [8].  A suitable  completion  of  these 

spaces  into  a Hilbert  space  and  construction  of  the  reproducing  kernel  was 
recently  accomplished  [1].  Since  these  spaces  are  made  up  of  functions  whose 
partial  derivatives,  up  to  a certain  order  in  each  variable,  are  absolutely 
continuous,  these  spaces  contain  spline  functions  in  two  variables,  and  the 
optimal  approximations  are  splines.  The  Sard  corner  spaces  have  the  property 
that  the  representers  reduce  to  products  of  functions  of  one  variable,  thus 
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simplifying  matters  somewhat.  Associated  with  the  Sard  space  Bpp  is 
a base  point  (a,b)  in  the  region  of  interest.  While  this  point  can  theo- 
retically be  anywhere,  practically  it  is  desirable  for  it  to  be  at  one  of 
the  corners  of  the  region  of  interest  (assumed  rectangular).  There  are  two 
reasons  for  this;  (i)  Computation  of  the  representers  is  simplified;  and 
( i i ) the  representers  have  continuous  partial  derivatives  in  x only  up 
through  order  p - 1 at  x = a and  in  y only  up  through  order  q - 1 at 
y = b , while  at  other  points  partial  derivatives  in  x and  y are  contin- 
uous up  through  order  2p  - 2 and  2q  - 2 , respectively.  The  second  reason 
is  the  primary  reason  we  assume  that  (a,b)  is  the  origin  and  that  we  are 
only  interested  in  (x,y)  points  in  the  first  quadrant. 

In  connection  with  the  previous  paragraph,  we  note  that  approximations 
in  Bpp  q-|  are  not  invariant  with  respect  to  translation  (unless  the  base 
point  (a,b)  is  also  translated),  nor  with  respect  to  stretching  or  shrink- 
ing of  the  coordinate  system.  The  author  had  previously  commented  that  they 
were  invariant  [3].  The  base  point  (a,b)  is  the  point  at  which  Taylor 
series  (with  remainder)  for  the  functions  exist,  and  the  approximations  are 
clearly  dependent  on  that  point. 

In  the  general  case,  the  reproducing  kernel  function  for  B^p  -|  is  of 
the  form  K(a,b;u,v,x,y)  = gp(a;u,x)gq(b;v,y)  where  for  a s u,x  , 

9_(a;u,x)  = (-l)p(x  - u)j2p_1)  + l l(u  - a)^(x  - a)^  + 
w i <p 

(-1 ) 1 (x  - a)(i+1)(u  - a)(P’1“1)}  . 

Here  the  notation  w^  means  wVi!  while 

!w  , w s 0 
0 , w < 0 . 
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For  the  case  we  consider,  with  a = b = 0,  p = q , and  where  the  func- 
tionals are  point  evaluations,  say  at  the  point  (x^.y^  , the  representers 
are  K^x.y)  = g^Ojx^  ,x)g^(0;yi.  ,y)  . As  in  what  follows  we  have 
here  used  simplified  notation  for  the  representer  associated  with  (x.,y..)  . 
We  will  consider  q = 1 and  q = 2 and  thus  note  that  g^(0;x..;x)  = 1 + x - 
(x  - xi)+  and  g2(0;xi;x)  = 1 + x.x  + ^ x.x2  - x3  + ^x  - x. . 

One  observation  about  optimal  approximations  is  now  in  order,  and  should 
lead  to  increased  interest  in  their  use  for  approximations  where  the  data  is 
irregularly  spaced.  The  data  points  themselves  generate  the  r^nresenters , 
and  hence  a set  of  basis  functions  for  the  approximations.  One  does  not 
need  to  be  concerned  about  whether  or  not  the  basis  functions  have  the  inter- 
polation property  on  the  set.  Unlike  more  common  basis  functions,  e.g., 
polynomials,  the  representers  naturally  form  a linearly  independent  set  over 
the  data  points.  This  is  not  meant  to  imply  that  the  representers  of  point 
evaluation  functionals  are  well  suited  to  computation,  however.  We  treat 
this  problem  in  more  detail  in  later  sections. 

1.1  The  interpolation  problem 

The  underlying  problem  we  shall  be  considering  is  that  of  function 
approximation  by  interpolation  for  functions  of  two  or  more  variables.  The 
case  of  more  than  two  variables  is  a straightforward,  if  tedious,  general- 
ization, and  the  discussion  is  limited  to  two  independent  variables.  Assume 
that  the  points  (xk ,yk*zk^ * k = 1**-*»N  are  given.  No  assumptions  are 
generally  made  as  to  the  spacing  of  the  points,  although  in  some  instances 
we  will  consider  special  cases.  We  assume  that  if  i * k,  (x..,y..)  * (x^,y^)  • 

Approximation  problems  other  than  interpolation  are  treated  in  identical 
fashion.  One  obtains  the  same  coefficient  matrix  for  approximate  integration 
or  differentiation,  for  example. 
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2.0  Optimal  approximation  in  ^ 

In  this  section  we  will  consider  in  some  detail  the  problem  of  computing 
optimal  approximations  in  the  Sard  space  ^ . Again  we  emphasize  that 
the  region  of  interest  is  assumed  to  lie  in  the  first  quadrant  and  that  the 
base  point  (a,b)  is  taken  to  be  the  origin. 

2.1  Representers  of  point  evaluation  functionals  as  a basis 

As  noted  in  the  introduction,  the  representers  of  point  evaluation 
functionals  have  the  form  K.(x,y)  = [1  + x - (x  - x.).][l  + y - (y  - y-)+]  . 
These  functions  K.  are  continuous,  with  first  partial  derivatives  in  x 

J 

and  y which  have  jumps  at  x = x-  and  y = y.  , respectively.  In  each 

J J 

of  the  rectangles  [0,Xj]x[0,yj] , [O.XjMyj ,») , [x^  .“OxEO.yj]  , and 

[x. ,»)x[y . ,o°)  the  function  is  bilinear.  In  the  latter  rectangle  it  takes 
J J 

on  the  constant  value  (1  + x.)(l  + y.)  . 

J J 

Use  of  the  K.(x,y)  as  basis  functions  appears  to  be  suspect.  The 

J 

coefficient  matrix  (the  Gram  matrix)  is  symmetric,  but  casual  observation 

would  lead  one  to  suspect  it  is  not  particularly  well  conditioned.  We  shall 

see  that  it  is  better  than  the  author's  inclination  toward  it.  The  system 

of  equations  has  the  form 
N 

(1)  E A.K  (x.,y.)  = z.,  i = 1 N . 

j=l  J J 1 1 1 

Some  numerical  experiments  were  conducted  to  compute  the  condition  number 
(with  respect  to  max  row  sum  norm)  of  some  Gram  matrices  of  various  sizes. 

Points  were  generated  by  a random  number  generator  in  the  square  [0,10]x[0,10]  , 
the  Gram  matrix  was  formed,  and  the  condition  number  computed.  In  case  (i)  the 
points  were  allowed  to  be  anywhere  in  the  square.  In  case  (ii)  the  square  was 
subdivided  into  [/TT]  squares,  then  one  point  was  generated  at  random  in  each 
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smaller  square,  and  any  remaining  points  were  generated  at  random  throughout 
[0,10]x[0,10]  . In  cases  (iii)  and  (iv)  the  points  in  (i)  and  (ii)  were 
translated  to  the  square  [90,100]x[90,100]  . In  cases  (v)  and  (vi)  the 
points  in  (i)  and  (ii)  we  stretched  over  the  interval  [0,100]x[0,100]  . 

In  case  (vii)  points  with  integer  coordinates  in  [0,lG]x[0,10]  were  selected 
by  a random  number  generator.  These  cases  give  some  indication  of  condition 
numbers  and  the  effects  of  translation  and  stretching.  No  claim  is  made  that 
the  cited  results  are  representative  nor  that  an  exhaustive  set  of  calculations 
were  made.  Other  computations,  not  tabulated  here,  follow  the  same  trend, 
however.  The  results  of  the  calculations  are  tabulated  in  Table  1.  Mot  all 


combinations  were 
Case  N -*■ 

computed. 

10 

25 

50 

100 

(i) 

170 

3400 

6700 

67000 

(ii) 

420 

1500 

5700 

36000 

(iii) 

24UU 

48UUU 

1 dUUUU 

(iv) 

3600 

37000 

200000 

(v) 

190 

3200 

8600 

(vi) 

600 

6100 

6200 

(vii) 

2900 

15000 

Table  1:  Condition  Numbers  of  (Kj(x^,y^))  for  Bp  ^ 

One  can  make  several  observations  from  the  table.  First,  for  moderate 
values  of  N , even  up  to  50  or  100,  satisfactory  results  can  easily  be  obtained 
by  computing  with  the  representers  as  basis  functions.  Depending  on  the 
accuracy  required  in  the  computed  answer  and  the  precision  one  can  (or  is 
willing)  to  use  in  the  computations,  N could  be  quite  a large  number,  per- 
haps as  large  as  is  feasible  to  consider  for  a global  approximation.  Second, 
the  effect  of  stretching  on  the  condition  number  is  rather  mild.  Third,  the 
translation  of  the  points  away  from  the  base  point  increases  the  condition 
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number,  indicating  that  one  should  select  the  base  point  (a,b)  close  to  a 
corner  of  the  region  of  interest,  as  we  have  done. 

2.2  A set  of  basis  functions  with  compact  support 

It  is  possible  to  construct  a set  of  basis  functions  with  compact  support 
which  leads  to  a block  tridiagonal  system  of  equations,  where  each  block  is 
tridiagonal,  for  approximations  in  j-j  . Unfortunately  the  scheme  is 

practical  only  in  the  special  instance  that  the  data  points  lie  on  a grid. 

Such  a grid  exists  for  any  set  of  points,  but  ordinarily  only  one  point  lies 
on  each  horizontal  and  each  vertical  grid  line.  For  the  scheme  to  be  prac- 
tical it  is  not  necessary  that  each  grid  point  be  a data  point,  but  many 
points  should  be  on  most  grid  lines.  This  will  be  made  more  explicit  after 
the  development. 

We  must  alter  our  usual  notation  slightly.  Suppose  we  have  grid  points 

(x.j,y.),  i = 1 . . , n,  j = 1 m . We  assume  that  the  x.  and  y,  are  in 

increasing  order.  Denote  the  set  of  subscript  pairs  corresponding  to  data 

points  by  I . Then  corresponding  to  each  (k,2.)el  there  is  a known  function 

value  zkA  . Nondata  points  on  the  grid  will  be  denoted  as  (i,j)^I  , where 

we  will  always  assume  that  1 s i s n and  1 s j < m . 

It  is  easy  to  obtain  functions  of  one  variable  of  the  appropriate  form 

which  have  compact  support.  They  are 

n 

G. (x)  - £ (0;x  »x) , i * 1 ,2 n , 

i r=l  ir  1 r 

where  as  before  g1  (0;xr,x)  = 1 + x - (x  - xp)+  , 

1 + x2  1 

with  an  = a + x'^u---")  * ai2 = ~-”x7 

1 xi+!  ~ xi-l  1 

“ii-!  * xi_1  - x.  ’ aii  ' (x.  - x.+1)(x.“i  - x^  ’ aii+l  “ x.  - x.+] 
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1 


1 


*nn-r  xn_1  - xn  * ann  “ xn  - x^ 

and  a^r  = 0 if  |i  - r|  > 1 . For  convenience,  let  xQ  = 0 and  x +.j  = ® , 
then  we  note  that  G. (x.)  = 6..  as  well  as  G.(x)  * 0 only  over  the  interval 

I J i J I 

(x-j_i»  x-j+-|)  • Also  construct  the  dual  functions  Hj(y)  > with 

Hj(y)  = z Sjs  g-j  (0;ys>y)  . Note  that  the  and  are  linear  B-splines. 

To  satisfy  the  interpolation  requirements  in  terms  of  the  local  basis 

n m 

functions  G. (x)H .(y) ,we  obtain  the  equations  z z a.  ,G. (x, )H . (y  ) = 

J -j  =]  j = l 'J  1 K J * 

ak,z  = zk,z  » (k,z)  £ 1 • For  (M)  t 1 * the  products  g-,  (0,x|<,x)g1  (0;yz,y) 
cannot  appear  in  the  approximation.  Substituting  for  G.(x)  and  H ^ (y ) the 

n m 

approximation  becomes  z ' Z a..  G.(x)H.(y)  = z a.,  z a.  g,(0;x  ,x)  z 

i=l  J-l  1J  1 J i,j  1J  r lr  1 r s 

Sjs  ■ I «,J  £ «,r  8js  g1(0;xr.x)g,(0;ys ,y)  . 

I IjS 

We  then  set  the  coefficient  of  g1  (0;xk>x)g-|  (Oiy^.y)  for  (k.2.)  i I , 

equal  to  zero,  obtaining  as  the  system  of  mn  equations  for  the  a..  , 

* 0 

( ak,t  = zk,£  * € 1 

(2) 

( aij  aik  8j*  = o.  <k’*)  * 1 • 

If  we  order  the  equations  and  variables  in  some  logical  fashion,  say 
(1.1)....  (l.m),  (2,1),...,  (2,m),...,  (n,m),  for  m < n , the  fact  that 
«i  k = 0 it  1 1 - k | > 1 and  3^  * 0 if  |j  - t|  > 1 shows  that  the 
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resulting  coefficient  matrix  of  the  system  of  equations  is  of  the  form 


/ Tn 


T = 


'12 


T21  T22 

0 \ 

T23 

0 

^"nn-1 

T 

l 

/ 611 

e12 

0 

' a21 

622 

623 

Tij  °1j 

O O 

X 

X 

0 \ 

0 

I 

nn  / 


0 \ 

0 


rnm-1 


mm 


/ 


except  that  a row  of  T corresponding  to  (k,i)  e I is  replaced  by  a new 
row  with  zeros  everywhere  except  for  the  unit  diagonal  element. 

We  can  note  that  6 j j _ i and  Sjj+-|  are  nonpositive  and  except  for 
j = 1 and  j - m,  Sjj  = - Sjj_,  - , while  for  J = 1.  s„  > - s,2 

and  for  j = m , . Thus  each  block  T.j  is  diagonally 

dominant.  An  interesting  aside  is  that  the  system  (2)  looks  very  much  like 
the  system  of  equations  obtained  when  solving  Laplace's  equation  by  finite 
differenceson  a rectangular  grid,  where  a 9 point  approximation  to  the  Laplacian 
is  used.  Here,  of  course,  the  "boundary  conditions"  may  be  scattered  through- 
out the  region.  In  this  respect  the  coefficient  matrix  is  well  suited  to  the  use 
of  iterative  methods  for  solution  of  (2). 

2 

For  the  usual  case,  the  above  leads  to  a system  of  N equations,  and  the 
amount  of  computation  will  be  too  large  for  even  moderate  N.  Use  of  a block 
elimination  equation  solver  reduces  the  problem  to  one  of  repeated  solution 
of  n systems  of  equations  of  m equations  (usually  full,  but  not  always). 
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3 

Although  the  original  problem  can  be  solved  in  about  N /6  operations,  the 

3 4 

new  basis  requires  about  n m , or  about  N operations  if  m = n = N . 

One  should  point  out,  however,  that  after  obtaining  the  a.,  one  has  sufficient 

* J 

information  to  obtain  the  value  of  the  interpolating  function  by  bilinear  inter- 
polation since  a,,  represents  the  function  value  at  (x.,y.)  . 

• J ' J 

We  will  consider  the  case  in  which  there  are  several  data  points  on  each 

grid  line.  The  amount  of  work  required  to  solve  (2)  decreases  somewhat  as 

the  fraction  of  known  grid  values  increases,  however,  the  major  effect  is  that 

the  relationship  between  (n,m)  and  N changes.  Assume  that  the  fraction 

of  grid  points  at  which  data  is  known  is  p,  0 < p s 1 . The  total  number 

of  data  points  is  then  N = pmn  , and  the  number  of  operations  required  to 

3 3 3 

solve  (1)  is  about  (pmn)  /6  . For  (pmn)  /6  =s  nm  we  see  that 

p »6/n^  , or  p % ^ n~6/  . Representative  values  of  these  fractions  are 

given  in  Table  2,  along  with  the  total  number  of  equations  in  (1)  and  (2). 

If  p is  larger  than  the  listed  value  in  Table  2,  (2)  can  be  solved  in  fewer 

operations  than  (1). 


n 

P 

N 

mn 

4 

.721 

2.88m 

4m 

8 

.454 

3.63m 

8m 

15 

.299 

4.48m 

15m 

50 

.134 

6.69m 

50m 

100 

.0843 

8.43m 

100m 

1000 

.0182 

18.2m 

1000m 

Table  2:  Fraction  p of  grid  points  to  be 
known  for  comparable  number  of  operations  in 
solving  (1)  and  (2). 
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2.3  Basis  functions  which  are  zero  in  some  regions 

Because  the  representers  of  point  evaluation  functionals  are  constant  for 
large  enough  values  of  the  independent  variables,  a certain  linear  combination 
of  any  two  can  be  made  zero  for  large  enough  values  of  the  independent  vari- 
ables. In  particular,  consider  K^x.y)  and  Kj(x,y)  . The  function 
(1  + x.  )(1  + y^KjU.y)  - (1  + xJ.)(l  + yj)K.(x,y)  is  zero  for  x 2 max(xi,xJ.) 
and  y 2 max(y. ,y.)  . Thus  a new  set  of  basis  functions  with  zero  values  over 

* J 

part  of  the  region  of  interest  can  easily  be  constructed.  It  is  desirable  to 
first  order  the  data  points  in  terms  of  their  "distance''  from  the  origin.  It 
seems  reasonable  to  order  the  (x^^y^)  in  terms  of  nondecreasing  values  of 
(1  + *k)(l  + y^)  • This  is  not  the  only  ordering  which  can  be  used,  but  it 
carries  the  assurance  that  the  new  jth  basis  function  will  be  nonzero  at  the 
point,  and  also  has  an  added  benefit  we  will  discuss  later. 

Assume  that  the  data  points  are  ordered  so  that  pk  = (1  + xk)(l  + , 

k * 1.....N  is  a nondecreasinq  sequence. 


Then  define  L,(x,y) 

J 


P1Kw(x,y)  - pjtlK|(«.y) 


and 


LN(x,y)  = 


KN(x,y) 

lV*N,yN) 


= 1 


,N  - 1 


Then  the  Lj(x,y)  satisfy  L ^ (x^  ) = 1,  j = 1,...,N  and 

It  ★ 

L j ( x ,y ) = 0 for  x 2 x^  = max(Xj ,Xj+1 ) and  y 2 y^  = max(yj ,yj+1 ) , j = 1... 

Using  these  basis  functions,  the  interpolation  problem  requires  the  so- 
lution of  the  system 
N 

(3)  ^ AjLj(Xj,y^)3z^»i-l,...,N. 

0 
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Because  of  the  above  ordering  on  the  data  points,  there  is  some  possibility 
of  the  entry  Lj(x.,yi)  being  zero  for  j < i , that  is  an  element  below 
the  diagonal  of  the  coefficient  matrix. 

The  basis  functions  L-(x,y)  have  a property  which  is  desirable  and 

J 

which  arises  out  of  our  ordering  of  the  data  points. 

Proposition:  In  the  first  quadrant,  |L.(x,y)|  s 1 . 

J 

We  note  that  it  is  also  true  that  if  Pj+-|  < Pj  . then  |Lj(x,y)|  < 1 
except  along  one  of  the  rays  which  start  at  (x-,y.)  and  extend  horizontally 

J J 

to  the  right  and  vertically  upward.  The  proof  is  simple  and  will  not  be 
given.  It  consists  of  considering  the  value  of  L-(x,y)  at  all  points  where 

J 

the  first  partial  derivatives  are  discontinuous,  since  any  extrema  must  occur 
at  such  a point.  These  points  are  (0,0),  (x^.O),  (xj+1,0),  (0^),  (0,yj+1), 
(Xj.yj),  (xj»xj+i)»  (Xj+i’V  , and  ( Xj+i  ,yj+i ) . Typical  function  values 
are  shown  in  figure  1.  Recall  that  in  each  rectangle  L.(x,y)  is  bilinear  and 

J 

hence  determined  by  its  values  at  the  corners. 
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Certain  behavior  can  be  classified  further.  For  example:  (i)  if 


*j+1  > xj.yj+i  > yj  » then  |Lj(x,y)|  < 1 except  at  (x^)  ; (ii) 

if  Pj+i  = Pj  » then  Lj(x,y)  = 0 for  x ^ min(Xj ,Xj+1 ) and 

y < min(y  .,y  .+1)  and  Lj(xJ+1  .y^  ) = -1  . 

The  coefficient  matrix  of  the  system  (3)  has  its  largest  element  (in 

magnitude)  on  the  diagonal,  and  some  zeros  may  occur  below  the  diagonal. 

The  data  points  used  to  generate  Table  1 were  used  to  test  the  effect- 
iveness of  the  introduction  of  zeros  and  to  determine  the  condition  numbers 


associated  with  the  new  basis  functions.  The  results  are  shown  in  Table  3, 
with  the  number  in  parenthesis  indicating  the  number  of  leading  zeros  in 
the  matrix.  By  reordering  the  columns  of  the  matrix  it  is  sometimes  possible 
to  introduce  many  more  leading  zeros,  and  while  a scheme  of  this  sort  has  not 
been  implemented,  in  many  cases  it  would  substantially  reduce  the  number  of 
operations  required  for  solution  of  the  system  (3). 

The  system  (3)  does  not  have  a symmetric  coefficient  matrix,  and  unless 
approximately  30%  leading  zeros  are  introduced,  it  will  take  fewer  operations 
to  solve  (1)  than  to  solve  (3).  However,  the  condition  number  of  the  new 
coefficient  matrix  has  been  smaller  in  every  case  examined,  sometimes  by  a 
factor  of  15  or  more,  but  more  commonly  by  a factor  of  2-5.  Reordering 
columns  for  a maximum  number  of  leading  zeros  will  often  result  in  30%  or 
more  leading  zeros,  based  on  some  hand  computations  where  no  effort  was  made 
to  obtain  the  maximum  number  of  leading  zeros. 
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Case  N -*• 

10 

25 

50 

100 

(i) 

36(4) 

250(20) 

1800(261) 

18000(535) 

(ID 

32(18) 

460(32) 

2000(162) 

12000(582) 

(ill) 

1800(6) 

2600(39) 

30000(315) 

(iv) 

160(18) 

3300(62) 

19000(288) 

(v) 

23(2) 

150(20) 

1400(195) 

(vi) 

28(18) 

390(16) 

1500(105) 

(vii) 

620(21) 

6200(146) 

Table  3:  Condition  Numbers  of  (Lj(xi,yi))  and  Number  of 

Leading  Zeros  (in  Parenthesis) 


3.0  Optimal  Approximation  in  Br-  .-i  . 


The  problem  of  computing  optimal  approximations  in  the  Sard  space 
B 2]  1S  somewhat  more  difficult  than  for  ^ . The  reproducing  kernel 
functions  are  seen  to  be  piecewise  bicubic  functions,  reducing  to  bilinear 
functions  for  sufficiently  large  values  of  the  independent  variables.  We 
shall  investigate  the  feasibility  of  extending  the  results  of  the  previous 
section  to  2]  in  secti°n* 


3.1  Representers  of  point  evaluation  functionals  as  a basis 

The  representer  of  the  point  evaluation  functional  at  the  point  (x.,y.)  is 

2 3 J J 

Kjlx.y)  ■ [’  ♦ y + * - f *3  * 5<*  - ♦ yjy  * jy/  y - jy3  - j<y  - yj)' 


This  function  is  cubic  in  x for  0 $ x s x . , and  linear  in  x for 
x > Xj  , and  the  dual  holds  in  y . These  functions  Increase  rapidly  since  the 
point  x * Xj  is  the  inflection  point  of  the  cubic  in  x , and  thus  when  the 
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function  is  linear  in  x , it  has  slope  the  same  as  the  maximum  slope  of  the 
cubic.  Because  of  this,  the  Gram  matrix  is  not  well  conditioned.  Following 
a similar  path  to  that  taken  in  the  previous  section,  some  condition  numbers 
for  the  Gram  matrix  were  computed  for  some  sets  of  randomly  generated  points. 
The  results  are  given  in  Table  4.  The  point  description  column  refers  to  the 


descriptions  in 

section  2.1 

• 

Case  N -*■ 

10 

25 

50 

100 

(i) 

3.94-104 

3.51-1 07 

8.27-107 

4.58*10 

(11) 

3 . 89  - 1 05 

1 . 89-1 06 

4 . 47  * 1 09 

9.95*10' 

Table  4: 

Condition  Numbers  of 

(Kj(xi,y.)  for 

Bf2,2l 

The  observations  which  we  wish  to  make  are  that:  (1)  One  will  quickly  be 
in  numerical  trouble  in  Real  *4  on  the  IBM  360,  and  (2)  While  the  condition 
number  is  large,  meaningful  computations  can  be  done  in  Real  *8. 


3.2  A set  of  basis  functions  with  compact  support 

A similar  construction  for  2-j  as  was  pursued  for  Bp  ^ in  section 

2.2  can  be  done.  There  is  some  question  as  to  what  conditions  should  be  imposed 

on  Gr  G2,  G3  , G4,  Gn_3,  G|N_2,  Gn_j,  and  GN  , (as  well  as  the  corresponding 

dual  functions  H. ),  but  several  reasonable  options  are  open.  In  the  general 
case  one  wants  V„  ’ j,1r  »«  W * ' • «.«  - ° • 


x i xi-2  or  x - xi+2  ' proceedin9  1n  similar  fashion  one  will  obtain  a system 
of  equations  which  in  block  form  has  5 non-zero  blocks,  each  block  being  a 
square  matrix  with  5 non-zero  elements  per  row.  Numerically  this  is  somewhat 
more  complicated  than  before,  but  there  are  instances  where  it  could  be  useful. 
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No  details  of  the  construction  have  been  carried  out  here,  but  the  general- 
ization is  straightforward.  Also  note  that  the  G.  and  here  are  cubic 
B-splines. 


3.3  Basis  functions  which  are  zero  in  some  regions 

It  seems  natural  to  be  able  to  extend  this  idea  to  Bp2  2]  ’ Tbe  bas1s 
functions  Kj(x,y)  are  bilinear  for  x > and  y >_ y ^ , thus  It  seems 

possible  a certain  linear  combination  of  them  could  be  made  identically  zero 
to  the  right  and  above  all  points  (xj ) associated  with  those  five  basis 
functions.  Proceeding  in  the  obvious  fashion,  ordering  the  points  (x^.y^) 
by  some  rule,  we  then  wish,  for  j < N - 4 , to  construct  functions 


Lj(x,y) 


j+4 

I.  Yjk  Kk(x,y)  such  that 
k-J 


0 


for  (xJ+4,  xj+3.  xJ+2,  xj+, 

and  y xyj.max  (yJ+4,  ,J+3.  yj+2.  yj+] 


and  LjfXj.yj)  = 1 . Unfortunately,  if  the  points  (xj+4,yj+4) , (xj+3,yj+3), 
(xj+2»yj+2)»  (Xj+^.y3+1)  lie  on  any  bilinear  curve,  this  system  of  equations 
generally  has  no  solution.  Thus  the  ordering  imposed  earlier  would  at  least 
have  to  restrict  one  away  from  four  successive  points  lying  on  a bilinear 
curve.  In  general,  this  is  not  possible. 

Because  of  the  rather  more  restrictive  region  where  the  function  is  zero, 
less  benefit  is  likely  to  accrue  anyway.  In  addition,  even  when  one  can  con- 
struct such  sets  of  basis  functions,  it  is  not  possible  to  bound  the  function 
Lj(x,y)  as  in  the  B^  ^ case,  and  in  particular  it  cannot  be  bounded  by  one. 
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4.0  Conclusions 


— ^ 

This  investigation  determined  that  use  of  the  representers  of  point 
evaluation  functionals  can  be  used  as  a basis  for  some  problems  without  en- 
countering severe  computational  problems.  For  smooth  approximations  this  is 
probably  not  generally  true,  however.  In  addition, ^because  the  computational 
burden  for  global  approximations  is  likely  to  be  quite  large,  it  is  the 
author's  opinion  that  local  approximations  must  be  investigated  for  smooth 
interpolation. ^.The  time  is  perhaps  propitious  for  an  investigation  into  the 
underlying  mathematical  basis  for  some  previously  suggested  schemes  for  local 
smooth  interpolation. 

The  use  of  optimal  approximations  in  Sard  corner  spaces  for  the  inter- 
polation of  irregularly  spaced  data  results  in  an  approximation  which  has 
discontinuities  along  the  lines  parallel  to  the  axes  through  each  data  point. 
This  would  seem  to  the  author  to  be  unnecessarily  complicated,  and  the  author 
intends  to  investigate  global  approximations  in  which  the  discontinuities 
are  less  numerous  and  local  in  character  rather  than  extending  alona  lines  to 
infinity. 
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